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1.  Summary 


The  goal  of  Phase  1  of  this  research  was  to  complete  development  and  verification  of,  as  a 
proof  of  concept,  an  adaptive  version  of  the  community  Fifth-Generation  NCAR/Penn  State 
Mesoscale  Model  MM5  including  a  coupled  advanced  hybrid  Large  Eddy  Simulation/Reynolds- 
Averagcd  Navier-S tokos  (LES/RANS)  turbulence  model  for  direct  prediction  of  optical  scale 
turbulence.  Upon  successful  completion  of  Phase  I.  the  adaptive  and  turbulence  algorithms 
developed  and  installed  in  MM5  were  to  be  installed  in  the  Advanced  Research  version  of  the 
next- generation  mesocale  prediction  Weather  Research  and  Forecasting  model  (WRF-ARW) 
and  also  verified  using  standard  reference  cases.  This  report  describes  progress  toward  these 
goals  and  will  describe  the  technical  advances  and  accomplishments. 

In  order  to  resolve  atmospheric  features,  such  as  gravity  waves,  to  the  degree  necessary 
for  accurate  optical  turbulence  prediction  and  to  ensure  maximum  benefit  from  the  advanced 
hybrid  LES-RANS  turbulence  model,  the  Euler  equations  in  MM5  are  transformed  to  a  gen¬ 
eral  curvilinear  coordinate  system  and  the  NCSU  dynamic  solution  adaptive  grid  algorithm 
(DSAGA)  is  installed.  This  adaptation  algorithm  reduces  grid  spacing  locally  by  adjusting 
the  position  of  grid  nodes  without  increasing  the  total  number  of  nodes  in  the  grid.  This 
method  of  adaptation  is  usually  referred  to  as  r-  refinement.  The  node  relocation  is  in  re¬ 
sponse  to  a  user-defined  weight  function  distribution  calculated  from  user  selected  features 
of  the  solution  Grid  adaptation  is  available  in  the  entire  domain  when  nesting  is  not  used, 
in  the  entire  lowest  nest  or  in  a  subset  of  the  lowest  nest. 

A  four-equation  hybrid  LES/RANS  turbulence  model  has  been  developed  and  incorpo¬ 
rated  into  MM5  to  provide  a  new  approach  for  the  calculation  of  the  index  of  refraction 
structure  function  Cft,  a  quantitative  measure  of  atmospheric  optical  turbulence.  These  four 
equations  are  used  to  model  the  turbulence  kinetic  energy  (A-),  the  variance  of  vorticity  (C). 
the  variance  of  potential  temperature  ( 6 "2)  and  its  dissipation  rate  ( fo ),  respectively. 

The  improvement  in  prediction  of  optical  turbulence  by  incorporating  DSAGA  arid  the 
hybrid  LES/RANS  model  in  MM5  is  documented  by  comparison  with  the  models  of  Dewan 
and  Jackson  applied  to  unmodified  MM5  output.  By  comparison  with  two  sets  of  observa¬ 
tions.  the  numerical  results  suggest  that  the  adaptive  MM5  C/  predictions  approach  binned 
observation  variation  in  regions  where  the  resolution  has  been  improved  by  DSAGA.  with 
little  additional  computational  resources.  Selected  comparisons  are  included  in  the  results 
section. 

The  development  and  verification  of  the  adaptive  MM5  is  documented  in  three  proceed¬ 
ings  publications  (References  1-3)  and  in  an  article  accepted  for  publication  in  the  A1AA 
journal  (Reference  4).  Reviews  of  the  research  were  given  at  three  AIAA  Aerospace  Sci¬ 
ences  Meetings,  one  American  Meteorological  Society  meeting,  and  three  workshops,  two  at 
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Hanscom  AFB  and  one  at  Arizona  State  University. 

The  NCSU  personnel  who  participated  in  this  project  were  D.  Scott  McRae.  Professor, 
PI;  Hassan  A.  Hassan.  Professor,  Co-PI:  Xudong  Xiao,  Research  Assistant  Professor,  Code 
Developer;  and  Yih-Pin  Liew,  Post-Doctoral  Research  Associate,  code  developer. 

The  second  phase  of  the  project  to  install  DSAGA  and  the  hybrid  model  in  WRF-ARW 
was  not  complete  at  project  termination.  The  primary  code  developer,  Dr.  Xiao,  left  the 
project  just  as  work  began  on  the  WRF  version.  Due  to  few  current  doctoral  graduates  with 
high-level  CFD  code  development  skills,  compounded  by  visa  restrictions  for  some  otherwise 
qualified  candidates,  eleven  months  and  an  international  search  were  required  to  find  and 
employ  a  qualified  replacement  for  Dr.  Xiao.  Although  the  WRF  version  is  not  complete, 
development  is  being  continued  under  separate  DOD  funding  and  will  therefore  be  available 
for  operational  use  when  completed. 

2.  Introduction 

Atmospheric  optical  turbulence  refers  to  the  range  in  the  turbulence  spectra  in  the  atmo¬ 
sphere  that  causes  significant  fluctuations  in  the  refractive  index  of  air.  These  fluctuations 
affect  the  propagation  of  optical  signals  and  images  by  random  refraction  which  can  result 
in  a  reduction  of  optical  signal  effective  power,  range,  and  coherence,  and  can  also  degrade 
image  quality.  The  primary  sources  of  high  altitude  atmospheric  optical  turbulence  are  grav¬ 
ity  wave  breaking  and  jet  streams.  Gravity  waves  arise  from  a  number  of  sources  including 
topography,  convection,  and  wind  shear.  The  turbulence  layers  resulting  from  these  gravity 
waves  have  thicknesses  of  the  order  of  tens  of  meters  and  may  extend  for  many  kilometers 
in  the  horizontal  direction. 

As  noted  by  Jumper  and  Belaud  (Reference  5),  considerable  observational  effort  has 
been  expended  with  the  goal  of  expanding  available  databases  of  the  refractive  index  in 
the  atmosphere,  and  average  profiles  have  been  obtained  for  selected  locations.  However, 
measuring  directly  the  detailed  state  of  the  atmosphere  for  a  wide  range  of  conditions, 
locations  and  times  is  not  practical.  Thus,  there  is  a  need  to  develop  simulation  tools  needed 
to  predict  optical  turbulence  for  given  initial  and  boundary  conditions.  These  simulation 
tools  can  then  be  initialized  with  local  conditions  and  observations  and  used  to  provide 
tactically  useful  predictions  to  meet  DOD  requirements. 

Examination  of  optical  turbulence  predictions  obtained  by  post  processing  distribu¬ 
tion  from  the  output  of  current  mesoscale  numerical  weather  prediction  (NWP)  codes  reveals 
two  fundamental  problems.  First,  the  numerical  resolution  is  not  adequate  to  resolve  gravity 
waves  or  the  scales  of  optical  turbulence.  Second,  the  optical  turbulence  parameterization 
via  Tatarskii  (C%)  is  based  on  methods  developed  for  Large  Eddy  Simulation  (LES)  scale 
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resolution  and  thus  is  not  well  suited  for  the  much  coarser  mesoscale  model  resolution.  Both 
of  these  NWP  code  problems  are  addressed  in  the  current  work;  first,  by  installing  a  dy¬ 
namic  solution  adaptive  grid  algorithm  (DSAGA)  to  improve  resolution  of  the  atmospheric 
phenomena  leading  to  optical  turbulence  and  second,  by  extending  a  physically  derived  tur¬ 
bulence  closure  (the  NCSU  k-C  model)  to  provide  directly  the  structure  functions  required 
to  predict  optical  turbulence.  Both  of  these  techniques  were  developed  for  aerospace  appli¬ 
cations  under  DOD  and  NASA  funding.  The  modification  of  these  techniques  for  NWP  use 
is  a  technology  transfer  example  and  is  an  unanticipated  benefit  resulting  from  the  original 
research  funding.  Progress  in  the  application  and  further  development  of  these  advanced 
techniques  for  optical  turbulence  prediction  has  been  reported  in  contract  reports,  and  in 
proceedings  and  an  article  by  Xiao,  et  al  (References  1  4). 

The  index  of  refraction  structure  function,  C2,  is  the  most  relevant  quantity  for  optical 
propagation  predictions.  Current  practice  for  the  calculation  of  G2  employs  mesoscale  NWP 
codes,  such  as  MM5  (Reference  6).  These  codes  forecast  relevant  atmospheric  quantities  such 
as  wind  velocity,  temperature,  relative  humidity,  etc.,  on  typical  grid  spacings  in  the  range 
4  30  km  in  the  horizontal  direction  and  0.3  1.0  km  in  the  vertical  direction.  The  vertical 
spacing  is  further  refined  in  the  planetary  boundary  layer  (PBL).  As  indicated  earlier,  optical 
turbulence  scales  are  of  the  order  of  tens  of  meters.  Therefore,  a  parameterization  of  C*  is 
employed,  which  gives  C2  in  terms  of  flow  properties  output  from  the  specific  mesoscale 
weather  prediction  code  being  used. 

The  expression  for  C„  for  optical  wavelengths  near  the  visible  range'  is  given  by  (Refer¬ 
ences  7.8) 

0*  =  f  76  x  10-»T  V  Ci  (1) 

where  P  is  the  pressure,  in  Pascals,  T  and  0  are  the  temperature  and  potential  temperature, 
respectively,  in  degrees  Kelvin,  and 

Cl  =  a2e„A1/3.  a2  =  2.8  (2) 

where  e.y  is  the  dissipation  rate  of  the  variance  of  potential  temperature  and  f  is  the  dissipat ion 
rate  of  the  variance  of  velocity  or  turbulent  kinetic  energy  (TKE). 

The  parameterization  developed  by  Tatarskii  (Reference'  7)  serves  as  the  basis  for  current 
statistical  models.  According  to  his  theory.  Cl  can  be  expressed  as 

d=2.8irg)2  p) 
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where  L$  is  the  “outer  length  scale”  and  2  is  the  vertical  direction.  Current  statistical 
models  (References  9, 10)  attempt  to  parameterize  L0  in  terms  of  altitude  and/or  velocity 
and  potential  temperature  gradients.  On  the  other  hand,  Jackson’s  model  (Reference  11)  uses 
the  temperature  lapse  rate  ( dT/dz )  to  parameterize  Lq  because  it  provides  better  correlation 
with  thermosonde  data.  Other  models  were  developed  by  Bougeault  et  al.  (Reference  12) 
and  Tunick  (Reference  13)  for  use  in  the  PBL.  Currently,  a  common  procedure  is  to  use 
outputs  from  mesoscale  models  as  inputs  to  the  above  models  to  calculate  C%. 

In  order  to  avoid  the  set  of  simplifying  assumptions  that  led  to  Eq.  (3)  and  current  mod¬ 
els  (References  9-13),  the  Euler  equations  in  MM5  are  supplemented  by  the  TKE  and  an 
equation  for  its  dissipation  rate.  vQ  ,  where  v  is  the  kinematic  viscosity  and  Q  is  the  variance 
of  vorticity,  or  enstrophy  (Reference  14);  and  an  equation  for  the  variance  of  potential  tem¬ 
perature,  6 "2,  and  its  dissipation  rate,  e#  (Reference  15).  These  four  equations  were  derived 
directly  from  the  full  Navier-Stokes  equations  and  thus  retain  true  relevant  physics.  As  a 
result,  it  is  now  possible  to  calculate  C „  from  Eqs.  (1)  and  (2)  without  further  simplifying 
assumptions.  Further,  in  order  to  account  for  the  scales  that  are  not  resolved  by  LES,  a 
hybrid  Large  Eddy  Simulation/ Reynolds  Averaged  Navier-Stokes  (LES/RANS)  solver  (Ref¬ 
erence  16)  has  been  incorporated  into  MM5.  The  RANS  model  provides  the  desired  subgrid 
scale  model  needed  in  LES  calculations. 

In  order  to  resolve  atmospheric  features,  such  as  gravity  waves,  to  the  extent  necessary 
for  accurate  optical  turbulence  prediction  and  to  ensure  maximum  benefit  from  the  advanced 
hybrid  LES-RANS  turbulence  model,  the  Euler  equations  in  MM5  are  transformed  to  a  gen¬ 
eral  curvilinear  coordinate  system  and  the  NCSU  dynamic  solution  adaptive  grid  algorithm 
(DSAGA)  is  installed.  This  r-  refinement  algorithm  reduces  grid  spacing  locally  in  response 
to  a  user-defined  weight  function  distribution  calculated  from  selected  features  of  the  solution 
and  does  so  without  increasing  the  total  number  of  nodes  in  1  he  grid. 

The  remainder  of  this  report  gives  the  technical  details  of  the  implementation  of  DSAGA 
and  the  hybrid  model,  and  compares  example  predictions  with  observations  carried  out  at 
Holloman  Air  Force  Base.  Conclusions  are  drawn  and  recommendations  made  for  future 
work. 
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3.  Technical  Approach  and  Modifications  to  MM5 

3.1.  Governing  Equations  and  Transformations 


The  non-hydrostatic  governing  equations  in  MM5  (Reference  G),  defined  in  the  .r.  y, 
a  coordinate  system,  are  transformed  in  all  three  dimensions  to  a  uniform  computational 
coordinate  system,  using  the  chain-rule  according  to: 


r  =  t  (4) 

&  =  £i{x,y,a,t),  i  =  1,2.3 


where  a  is  the  nondimensional  pressure  coordinate.  The  resulting  equations  have  the  form 


au  d\j_d^  .  2dFd^  dG oz, 

dr  ^  0£i  dt  ^  d£j  dx  d£,  dy  da 

where 


(5) 


U  =  \])*p',p*u.p*v,  p*  w.  p*T]  J . 

E  =  ^\p*p',p*u.p*v,p*w,p*T]T,  (6) 

F  =  —  \ifp\p*u,p*v>p*w,p*T}7\ 

rn 

G«  r  *  /  *  *  *  *  7* 

—  a\p  p  ,  p  u, p  v.  p  wyp  /  j  , 


in  is  the  map  scale,  p*  the  reference  pressure,  and  p'  the  pressure  perturbation,  T  the  tem¬ 
perature,  and  u.  v  and  w  are  velocity  components  in  the  x-,  .(/-and  e-direction,  respectively. 
All  other  terms,  such  as  pressure  gradient,  Coriolis  force,  and  gravity  terms  are  included  in 
S;  see  Reference  6  for  more  details.  The  above  equations  are  discretized  in  the  Arakawa-B 
(Reference  17)  type  staggered  grid,  using  the  same  finite  difference  stencils  as  MM5.  e.g..  the 
stencils  used  in  the  ./’-direction  in  MM5  are  applied  to  the  G  direction  here.  These  equations 
are  also  solved  using  the  leap-frog  scheme.  In  order  to  obtain  accurate  discretization  in  the 
curvilinear  staggered  grid,  three  sets  of  metric  derivatives  are  calculated  to  be  consistent 
with  the  differencing  of  flow  variables.  The  variables  and  metric  derivatives  are  defined  at 
three  different  locations:  the  cell  center  for  p'  and  T,  the  center  of  cell  edges  for  u  and  v  and 
the  center  of  cell  surface  for  w  and  a. 

In  MM5,  in  order  to  remove  the  limitation  on  the  time  step  due  to  small  mesh  spacing  in 


the  vertical  direction,  the  following  two  coupled  equations  for  w  and  p'  are  solved  implicitly: 


dj]_ 

dt 


dw 

dt 


Po9  dp'  gjj_ 
pp*  da  7  p 


PoTW  dw 
p *  da 


-  pogw 


Su, 


S„ 


(7) 

(8) 


This  results  in  a  tridiagonal  system  along  the  a-direction  for  w  in  the  uniform  mesh,  which 
can  be  solved  directly.  In  the  computational  plane,  the  a  variation  is  transformed  as  dff  = 
& i<Td-^i,  where  stands  for  |j^.  Therefore.  Eqs.  (7)  and  (8)  must  be  solved  iteratively.  The 
iteration  scheme  chosen  is  as  follows: 


1.  p'W  =  pn ,  =  wl ; 

2.  solve  the  following  system  iteratively: 


u;(,+1)  —  wl 
A t 


PP*  d£ 3  7  P 


(9) 


p/(*+1)  _  Po£W  ,dw  (i+n 

At  p *  [d£ 3  <:i 


.  -  =  s,  +  «522|(^)«'>6,  +  (#)<%,] 


P* 


'  ^2 ' 


(10) 


3.  when  converged,  p,t+1  =  p'h+1)  and  u,<+1  =  it/ 


1  +  1) 


where  t  and  t  +  1  denotes  two  time  levels,  respectively,  “(0)”  the  initial  values  and  i  the 
number  of  iterations.  For  simplicity,  p'^  and  u/‘)  are  used  to  illustrate  the  algorithm. 
However,  they  are  implemented  as  an  averaged  value  of  p'{l)  and  pn  ,  and  it/*)  and  wl, 
respectively  (Reference  6).  Note  that  Eqs.  (9)  and  (10)  are  in  a  form  similar  to  Eqs.  (7) 
and  (8)  so  that  the  tridiagonal  solver  can  be  applied  to  the  ^-direction.  To  implement  this 
semi-implicit  solver  in  transformed  space,  an  outer  loop  is  added  to  the  original  MM5  loop 
to  update  the  right  hand  side  of  the  above  equations. 
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3.2.  The  Dynamic  Solution  Adaptive  Grid  Algorithm  (DSAGA) 

The  DSAGA  adaptive  grid  algorithm  as  currently  installed  in  the  MM5  is  a  modification 
of  the  NCSU  developed  algorithm  and  was  first  described  in  Reference  1.  In  order  to  install 
the  adaptive  grid  algorithm,  the  non-hydrostatic  governing  equations  in  MM5  defined  in  the 
x,  y.  (t  coordinate  system  are  transformed  to  a  general  curvilinear  coordinate  system.  The 
resulting  equations  were  discretized  in  the  Arakawa-B  (Reference  17)  staggered  grid  and 
solved  using  the  leajr-frog  scheme.  Since  the  flow  properties  are  stored  at  three  different 
locations  in  the  Arakawa  grid,  three  sets  of  grid  metric  derivatives  are  calculated  that  define 
the  transformation  to  the  curvilinear  grid  as  noted  in  the  preceding  section.  The  DSAGA 
algorithm  then  performs  /'-refinement  adaptation  in  order  to  increase  resolution  of  solution 
features  by  relocating  grid  nodes  in  a  solver-independent  manner.  A  weight  function,  based 
on  user-selected  solution  properties,  is  calculated  and  then  used  in  a  center  of  mass  scheme 
to  control  the  grid  node  relocation.  In  the  work  presented  here,  the  weight  function  is  based 
on  the  distribution  of  the  magnitude  of  vortieity: 

ir  =  |Vxi/|.  (ii) 

where  1  is  the  velocity  vector.  As  has  been  shown  in  Reference  1,  this  weight  function  is 
a  good  indicator  for  gravity  wave-breaking  and  also  is  significant  in  other  regions  of  large' 
shear.  The  weight  function  due  to  vortieity  associated  with  each  cell  is  treated  as  the  mass 
of  the  cell.  Then  the  cent  cr-of- mass  algorithm  determines  the  new  grid  by  clustering  nodes 
where  the  magnitude  of  vortieity  is  large,  thereby  increasing  resolution  in  regions  of  large 
shear  (e.f.  Reference  1)  for  details.  The  basic  grid  adaptation  procedure  using  DSAGA  can 
be  outlined  as  follows: 

1.  Calculate  the  weight  function  using  Eq(ll). 

2.  Restrict  the’  weight  function  by 

W  =  min(inax(ir,  uiU’ave).  02B  aw)  (12) 

where  irav(,  is  the  average  value  of  W  in  the  entire  domain  and  a  1  and  0*2  are  two 
coefficients  that  prescribe  the  floor  and  ceiling  values  of  the  weight  function,  respec¬ 
tively.  This  reduces  drastic  grid  node  movements  due  to  very  large  or  very  small  weight 
functions. 

3.  Smooth  the  weight  function  using  an  elliptic  smoother  to  reduce  the  “noise”  in  the 
weight  function  in  order  to  obtain  a  grid  with  smoother  distribution  of  cell  volume. 
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Then  rescale  the  resulting  weight  function  based  on  an  input  max-to-min  ratio  in  order 
to  bias  the  relocation  of  the  grid  toward  larger  or  smaller  movement-. 

4.  Obtain  the  new  grid  displacement  in  the  computational  domain,  (A£,  A r/,  Ac,),  using 
the  center-of-mass  algorithm. 

5.  Calculate  the  new  grid  node  locations  in  the  physical  coordinate  system,  and  redis¬ 
tribute  the  flow  solutions  to  the  new  locations  using  a  high  order  advection  scheme. 

By  examining  Eq(12),  it  is  clear  that  this  restriction  tends  to  adjust  the  maximum  and 
minimum  value  of  the  weight  function  according  to  the  average  of  all  weights,  i.e.,  if  the 
average  of  the  weight  functions  are  small,  then  both  ceiling  and  floor  value  in  Eq(12)  are 
small.  In  practice,  the  values  of  a\  and  «2  require  careful  tuning  on  a  case  by  case  basis. 
Therefore,  isolated  high  vorticity  regions  may  not  influence  the  function  sufficiently  using 
this  approach.  Since  the  purpose  of  this  “restriction”  is  to  remove  too  large  and  too  small 
values,  a  new  approach  based  on  the  n-th  largest  value  search  (Reference  18)  is  employed: 

W  =  min(max(W;  W„,%),  Wn2%),  (13) 

where  Wnx%  {W„2%)  is  a  weight  function  value  that  is  greater  than  n\%  (712%)  of  population 
and  less  or  equal  than  the  rest,  and  rq  <  112.  In  a  sequential  mode,  the  cost  of  this  searching 
algorithm,  in  terms  of  both  time  and  memory,  is  proportional  to  the  number  of  grid  nodes, 
so  it  is  as  efficient  as  the  approach  in  Eq(12).  Using  this  approach,  the  floor  and  the  ceiling 
values  are  not  dependent  on  the  average  value,  and  it  can  guarantee  that  there  are  (7? 2  — rq)% 
of  weight  functions  unaffected  by  this  restriction  procedure. 

In  the  solver  independent  version  of  DSAGA,  as  applied  here,  the  grid  relocation  process 
determines  the  grid  distribution  such  that  it  improves  resolution  of  the  solution  as  an  initial 
condition  for  the  next  solver  step.  However,  the  solution  on  which  this  relocation  is  based 
is  defined  only  on  the  previous  grid.  A  conservative  advection  scheme  is  then  used  to  redis¬ 
tribute  the  solution  to  the  new  grid  nodes.  This  step  is  crucial  for  preserving  the  accuracy 
of  the  solution  and  for  realizing  the  benefit  of  the  relocation. 

In  the  solution  redistribution  step,  an  advection  equation  is  solved  for  each  flow  variable. 
The  first  scheme  used  for  this  task  was  a  WENO  (Reference  19)  based  scheme  that  tended 
to  smooth  the  solution  and  thereby  did  not  provide  the  required  redistribution  accuracy.  A 
modified  Piecewise  Parabolic  Method(PPM)  scheme,  developed  by  Rider  et  al  (Reference 
20),  was  adopted  to  increase  redistribution  accuracy.  I11  the  PPM  framework,  the  single 
time  step  scheme  can  achieve  third  order  accuracy  in  both  the  temporal  and  the  spatial 
directions  due  to  integration  over  a  constructed  parabola  within  grid  cells.  Therefore  it 
is  more  efficient  than  the  multi-step  WENO  scheme.  In  Rider’s  approach,  the  base  PPM 
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(Reference  21)  is  combined  with  WENO  and  high  order  (non-monotone)  schemes.  A  median 
function  is  employed  to  select  the  optimal  cell  edge  values  for  constructing  the  parabola 
in  the  cell.  This  approach  maintains  uniformly  high  order  accuracy,  even  at  extrema.  A 
slight  change  was  made  on  the  high  order  component,  to  simplify  the  implementation,  while 
maintaining  similar  accuracy.  Instead  of  a  6  point  stencil  as  in  Reference  20.  the  following 
4-point,  stencil  is  used  for  the  edge  value: 

0J  +  1/2  =  +  4>j  + 1)  _  y}(<Pj+  2  +  0j- t)  (Id) 

where  the  0  is  the  variable  of  interest. 


3.3.  Turbulence  Closure 

3.3.1.  Review  of  Tatarskii’s  Theory 

In  order  to  contrast  the  proposed  approach  with  existing  approaches,  a  brief  review 
of  Tatarskii’s  theory  is  presented.  Tatarskii  assumed  that  turbulence  exists  in  a  state  of 
equilibrium  and  that  the  Richardson  number  is  zero.  The  first  assumption  requires  that 
production  and  dissipation  rates  of  TKE  and  0"2  are  equal,  while  the  second  ignores  the 
contribution  of  gravity  to  the  above  terms.  Thus,  the  above  assumptions  yield 


e  =  v, 


(15) 


and 


(lb) 


where  u,  and  «,  are  the  turbulent  eddy  viscosity  and  turbulent  diffusivity.  Using  Eqs.  (15) 
and  (lb).  Eq.  (2)  reduces  to 


V2 

0 


(-) 

L(ii)2J 

1/3 

'W\2 

dz) 


(17) 


1'lie  next  assumption  of  Tatarskii  pertains  to  ut  and  a,.  He  assumed  that  i/t  is  determined 
from  the  mixing  length  theory  of  Prandtl,  i.e.. 


-  ’Li 


chi 


dz 


=  Ll 


On 


dz 


(18) 


Moreover,  he  sets 


Vt 


1 


(19) 
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As  a  result  of  the  above  assumptions, 

Cl  =  2.SL J/3  (20) 

Existing  statistical  models  (References  9  11)  of  optical  turbulence  differ  by  the  manner 
in  which  Lq  is  modeled.  It  is  chosen  as  a  function  of  one  or  more  °f  §f ,  fz  an(^ 

3.3.2.  The  Hybrid  LES/RANS  Turbulence  Model 

In  this  work,  the  turbulence  is  modeled  using  a  hybrid  LES/RANS  model.  The  four 
equation  k-C,  model,  which  is  an  extension  of  the  standard  two-equation  model  (Reference 
14),  serves  as  the  RANS  component.  The  four  equations  account  for  the  turbulence  kinetic 
energy  (TKE)  or  the  velocity  variance,  A’,  the  variance  of  vorticity  or  enstropliy,  (\  which 
provides  the  dissipation  rate  (e  =  for  A-,  the  variance  of  the  potential  temperature,  0"2, 
and  its  dissipation  rate,  eg.  These  four  equations  were  derived  from  the  exact  Navier-Stokes 
equations  and  then  modeled  term  by  term  so  as  to  insure  that  relevant  physics  is  reflected 
in  the  resulting  model  equations.  The  LES  subgrid  scale  (SGS)  model  consists  of  the  k-A 
and  9"2- A  models,  where  the  A  is  the  grid  spacing  in  the  vertical  direction.  The  A:-A  model 
was  a  variant  of  the  one  equation  model  of  Yoshizawa  and  Iloriuti  (Reference  22),  which 
can  recover  a  Smagorinsky  model  in  LES  regions  when  production  balances  dissipation  in 
the  absence  of  buoyancy  effect.  The  current  O''2- A  takes  a  similar  form  as  the  k-A  model. 
Inclusion  of  equations  for  0"2  and  eg  makes  it  possible  to  calculate  the  turbulent  thermal 
diffusivity  as  part  of  the  solution.  Thus,  this  model  computes  the  turbulent  Prandtl  number 
directly. 

The  hybrid  scheme  results  from  the  combination  of  the  RANS  and  the  LES  components 
through  a  blending  function  (Reference  16),  T,  which  blends  the  turbulent  eddy  viscosity,  vt, 
turbulent  thermal  diffusivity,  at,  the  Adequation,  and  the  0"2-equation.  The  current  blending 
function  varies  from  0  to  1.  When  T  =  0,  the  hybrid  model  reduces  to  the  RANS  model, 
while  T  =  1  yields  the  LES  limit.  The  advantage  of  this  approach  lies  in  the  fact  that  all 
scales  are  taken  into  consideration.  The  scales  that  can  be  resolved  are  calculated  directly, 
while  the  ones  that  cannot  be  resolved  are  modeled.  If  the  resolution  is  insufficient  for  LES 
scales,  the  RANS  model  is  used. 

In  this  hybrid  scheme,  the  dissipation  rates  for  both  variances  consist  of  two  components: 
one  is  due  to  the  RANS  model,  and  the  other  is  due  to  the  LES  subgrid  model.  Therefore, 


in 


the  resulting  dissipation  rates  are  expressed  as 


(1-rK  +  r— 


a-3/'2 


(21) 


f  hybrid 


^ OJiybrid 


(22) 


where  vQ  is  the  dissipation  rate  for  A  in  the  RANS  component,  is  a  model  constant, 
and  Cp  is  the  speeihe  heat  at  constant  pressure.  With  these  two  dissipation  rates.  Cq  can 
be  directly  calculated  using  Eq.  (2),  instead  of  using  the  simplifying  assumptions  noted 
previously. 

In  order  to  take  the  buoyancy  effect  into  account,  the  heat  flux  0"u",  is  modeled  by 


(23) 


where  C)u,,  is  a  model  constant,  g  is  the  gravity  acceleration,  and  t(I  is  a  time  scale  deter¬ 
mined  by  0"'2  and  e o.hybrui ■  The  above  approach  is  based  on  a  similar  modeling  in  Garratt's 
monograph  (Reference  23). 

The  four  model  variables  are  initialized  by  assuming  that  the  initial  vt  equals  to  the 
Smagorinsky  eddy  viscosity,  Vtj.ES  —  Vl.hans  and  vt,Lt:s  ~  <*utANs-  The  model  equations 
are  discretized  and  solved  in  the  same  fashion  as  the  other  equations  in  MM 5.  Zero-gradient 
lateral  boundary  conditions  are  used  owing  to  the  absence  of  A\  (,.  B"‘2  and  e#  from  outer 
domains  in  the  standard  MM5.  The  model  equations  and  model  constants  used  are  given  in 
the  Appendix.  Details  of  the  derivation  and  model  constants  are  given  in  References  14  and 
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4.  Results  and  Discussion 


Selected  results  are  presented  in  this  section  to  illustrate  the  improved  forecast  pro¬ 
vided  by  the  NCSU  four  equation  hybrid  LES/RANS  turbulence  model  and  the  NCSU 
adaptive  algorithm  DSAGA  as  installed  in  the  modified  MM5.  More  complete  results  are 
available  in  Xiao,  et.al.  (Reference  1  4) 

Three  code  configurations  are  used: 

1.  standard,  unmodified  MM5  with  the  standard  MM5  turbulence  modeling. 

2.  The  adaptive  MM5  with  the  four  equation  turbulence  model  executed  in  fixed,  hori¬ 
zontally  uniform,  grid  mode  (referred  to  as  “uniform”  in  figures). 

3.  The  adaptive  MM5  with  the  four  equation  turbulence  model  in  full  three  dimensional 
adaptive  mode.  The  installed  capability  to  increase  the  number  of  levels  independently 
in  the  adaptive  domain  was  not  exercised  for  these  comparisons  as  tin1  initial  goal  was 
to  demonstrate  benefit  with  little  additional  computational  resources. 

Two  sets  of  thermosonde  measurements  were  used  during  the  project  to  provide  observa¬ 
tional  comparisons.  The  first  set  consists  of  three  thermosonde/radiosonde  profiles  collected 
in  Oct  2001  at  Vandenberg  Air  Force  Base  (Reference  24).  The  second  set  was  collected  in 
Oct  2003  at  Holloman  Air  Force  Base.  Results  for  all  three  code  configurations  will  be  com¬ 
pared  directly  with  the  Holloman  observation.  Comparisons  with  the  Vandenberg  sets  are 
included  in  Reference  2.  Balloon  number  and  launch  time  for  these  data  are  given  in  Table  1. 
The  high-resolution  observations  along  the  balloon  trajectory  were  binned  at  500m  vertical 
intervals  for  comparison  with  numerical  results.  In  addition  to  the  observational  data,  C„ 
profiles  obtained  from  Dewan’s  (Reference  10)  and  Jackson’s  models  (Reference  11)  using 
standard  MM5  output  (Reference  24)  and  using  adaptive  MM5  output  are  included  in  the 
comparisons.  For  these  two  models,  the  prediction  depends  solely  on  the  velocity  and  tem¬ 
perature  distributions  that  are  computed  by  the  standard  MM5  or  interactively  computed 
by  the  modified  MM5  using  the  hybrid  LES/RANS  model.  Since  the  LES/RANS  model 
interacts  with  the  solution,  the  Dewan  and  Jackson  results  for  C„  produced  by  these  two 
inputs  will  differ. 


Table  1.  Model  Initialization  time  and  Balloon  Launch  Time 


Model  Initialization  Time 

Balloon  # 

Balloon  Launch  Time 

Holloman  Case 

26-Oct-2003  1200  UTC 

1 

27-Oct-2001  0252  UTC 

However,  results  obtained  with  either  the  modified  or  original  MM5  cannot  be  expected 
to  reproduce  the  binned  thermosonde/radiosonde  data  even  if  nominal  vertical  spacing  is 
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Table  2.  Attributes  of  MM5  Grids  in  Holloman  case 


Nest 

1 

2 

Adaptive 

Horizontal  Grid  Spacing  (km) 

15 

5 

5  (initial) 

Horizontal  Grid  Size  (number  of  nodes) 

109  x  109 

163  x  163 

123  x  123 

Vertical  Grid  Size  (Levels) 

80 

similar  to  the  binning  range  in  the  stratosphere.  First,  numerical  dissipation  is  required  to 
maintain  stability  of  the  MM5  computation  and  is  also  used  to  damp  spurious  acoustic  modes, 
etc.  The-  collective  effect  of  these  dissipations  will  tend  to  reduce  the  amplitude  of  (or  even 
remove  entirely)  the  higher  frequencies  resolved  theoretically  by  the  mesh  spacing.  Second, 
the  absorptive  layer  condition  used  to  cancel  wave  reflections  at  the  upper  computational 
boundary  tends  to  damp  frequencies  in  the  upper  stratosphere.  These  effects  tend  to  smooth 
shear  layers  with  the  result  that  MM5  distributions  of  temperature  and  velocity  used  to 
obtain  C‘j  will  have  reduced  variation  and  structure  as  compared  with  observation.  This 
reduced  variation  will  also  influence  the  predictions  obtained  with  the  models  of  Jackson 
and  Dewan,  which  depend  on  these  distributions. 

Both  standard  and  adaptive  MM5  inputs  for  mesh  sizing,  nesting  and  initial  and  bound¬ 
ary  conditions  for  the  Holloman  case  were  those  used  previously  by  Ruggiero  (Reference  24) 
for  comparisons  with  observation  using  the  standard  MM5.  The  grid  attributes  for  standard 
and  adaptive  MM5  domains  are  given  in  Table  2.  The  adaptive  domain  is  nested  in  the  inner 
standard  MM5  domain  (Domain  2),  which  provides  lateral  boundary  conditions  needed  for 
the  adaptive  domain.  Initial  grid  resolution  in  the  adaptive  domain,  in  both  the  horizontal 
and  the  vertical  directions,  is  the  same  as  that  for  Domain  2  as  given  in  Table  2.  The  do¬ 
main  layout  showing  the  relative  location  of  the  grid  nest  and  adapted  region  is  shown  in 
Figure  1.  The  surface  normal  mesh  distribution  used  by  Ruggiero  (Reference  24)  results  in 
approximately  even  300m  vertical  spacing  in  the  stratosphere.  The  MM5  initialization  time 
for  this  case  is  1200  UTC  26-()ct-2003. 

During  computation,  grid  relocation  and  the  subsequent  solution  redistribution  are  con¬ 
ducted  every  three  Domain  2  time  steps.  The  time  step  in  the  adaptive  domain  is  calculated 
temporally  such  that  the  adaptive  solver  uses  the  maximum  time  step  that  the  CFL  condi¬ 
tion  allows  in  order  to  advance  to  the  same  time  level  as  in  the  outer  standard  MM5  domain. 
All  Cfx  resulting  from  the  adaptive  MM5  solver  are  calculated  using  the  four-equation  hybrid 
LES/RANS  model  described  in  Section  3.3. 

In  order  to  obtain  profiles  of  C*  and  other  quantities  along  the  balloon  trajectory,  three  di¬ 
mensional  flow  field  data  were  output  every  10  minutes  simulation  time,  and  TECPLOT™was 
used  to  interpolate  the  solutions,  spatially  and  temporally,  to  the  trajectory.  A  snapshot  of 
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two  vortical  grid  surfaces  from  the  three  dimensional  adaptive  grid  region,  with  the  weight 
function  contours  superimposed,  is  shown  in  Figure  2.  Note  that  none  of  the  adapted  grid 
surfaces  are  planar  since  the  adaptation  is  fully  three  dimensional.  The  balloon  trajectory 
is  projected  onto  this  snapshot  so  that  the  influence  of  the  wind  field  is  evident.  As  can 
be  observed  in  this  figure,  grid  nodes  are  clustered  in  regions  of  large  weight  function  (high 
vorticity).  In  regions  below  z=4  km  ,  the  vorticity  based  weight  function  was  set  to  the 
floor  value  (see  Eq.  13)  to  allow  concentration  of  adaptation  in  the  upper  atmosphere.  This 
view  could  also  have  been  superimposed  with  the  local  value  of  the  C%  prediction  in  order 
to  indicate  that  the  LES-RANS  turbulence  model  gives  this  prediction  everywhere  for  all 
model  run  times  as  a  direct  output. 

Figure  3  shows  a  comparison  of  the  potential  temperature  profiles,  where  good  agreement 
is  indicated.  This  figure  also  shows  that  the  tropopause  is  located  at  z  ~  11  km.  The  adaptive 
grid  algorithm  weight  function  is  biased  in  this  case  to  provide  more  adaptation  above  the 
tropopause,  where  accurate  prediction  of  is  required  for  this  application,  and  less  below. 
The  resulting  grid  is  refined  dynamically  in  the  region  between  z  ~  1 1  km  and  z  ~  19  km,  as 
indicated  by  grid  spacing  distribution  in  Figure  4.  An  earlier  grid  distribution  using  Eq.  12 
(instead  of  Eq.  13)  and  the  grid  distribution  used  both  in  the  unmodified,  fixed  grid  MM5 
and  as  input  for  the  adaptive  MM5  are  included. 

Figure  5  gives  a  comparison  of  C%  obtained  by  Dewan’s  and  Jacksons’  models  by  use  of 
output  from  the  unmodified,  fixed  grid  MM5  with  that  from  direct  output  of  Cj.  from  the 
current  adaptive  MM5  and  the  hybrid  LES/RANS  model.  As  a  result  of  grid  refinement, 
the  adaptive  MM5  C\  profile  shows  more  structure  in  the  upper  atmosphere  and  reproduces 
more  closely  the  large  C \  variations  shown  in  the  binned  observation  which  are  absent  from 
the  Dewan’s  and  Jackson’s  model  predictions  from  the  unmodified  MM5. 

Figure  6  compares  the  same  adaptive  LES/RANS  C\  profiles  shown  in  Figure  5  with  the 
adaptiveLES/RANS  code  in  fixed  grid  mode  (labeled  as  uniform).  These  results  demonstrate 
that  the  C/  prediction  provided  by  the  hybrid  LES/RANS  model  is  improved  in  both  high 
and  standard  resolution  grid  modes.  Figure  7  compares  the  meridional  wind  profiles  for 
results  shown  in  the  previous  Figures.  The  standard  MM5  profile  due  to  Ruggiero  (Reference 
24)  is  smooth,  indicating  that  most  strong  wind  shear  was  reduced  by  a  combination  of 
coarse  grid  spacing  and  dissipation  in  the  standard  MM5.  The  level  of  wind  shear  present 
in  the  binned  observation  is  more  closely  approximated  in  the  current  adaptive  results, 
indicating  that  the  overall  effective  dissipation  is  reduced  by  adaptation  and  the  four-equation 
turbulence  model.  The  misalignment  of  extrema  between  the  adaptive  solution  and  the 
observational  data  may  be  due  to  the  sponge  layer  on  the  upper  boundary,  which  may  cause 
a  boundary  layer  like  behavior  near  the  top  boundary  and  displace  the  extrema  locations 
towards  the  surface.  This  misalignment  of  extrema  is  also  apparent  in  the  C„  profiles. 
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Although  it  is  clear  from  the  above  noted  figures  that  the  displacement  of  the  results  by 
the  boundary  condition  implementations  preclude  any  meaningful  application  of  standard 
accuracy  measures,  examination  of  statistical  measures  may  indicate  progress  toward  the 
goal  of  direct  verification.  Table.  2  compares,  quantitatively,  the  performance  of  different 
schemes  in  both  low  and  upper  atmosphere  in  terms  of  the  mean  and  the  standard  deviation 
of  Cf.  The  mean  of  Cf  is  defined  as: 


and  the  standard  deviation  of  C f  is  defined  as 


where  Z\  and  Z2  are  the  lower  and  upper  bounds  of  altitude.  The  mean  describes  the 
integral  of  Cf  through  the  air,  which  is  an  useful  quantity  in  relevant  calculations  for  optical 
propagation  (Reference  9).  while  the  deviation  can  be  used  to  describe  the  variations  of  C‘j. 
Not e  that  the  Cf  presented  were  collected  along  the  balloon  trajectory.  However,  numerical 
experiments  show  that  110  significant  difference  exists  in  these  statistical  measures  between 
profiles  along  a  vertical  line  through  the  observation  site  and  along  the  actual  trajectory, 
since  the  balloon  trajectory  crosses  only  a  few  grid  cells  horizontally.  Therefore,  we  choose 
to  integrate  Cf  in  the  ^-direction. 

Table  3.  Comparison  of  mean  and  standard  deviation  of  at  lower  (4  km  <  ;  <  10  km)  and 
upper  (10  km<  z  <20  km)  atmosphere 


Models 

Lower  At 

mean 

mosphere 

<T 

Upper  At 

mean 

mosphere 

<7 

Observation 

5.09  x  10“18 

5.77  x  nr18 

1.99  x  HU17 

3.52  x  10" 17 

Dewan 

2.74  x  ltr17 

3.64  x  10"17 

7.61  x  HU18 

1.26  x  10"17 

Jackson 

3.51  x  10  18 

4.60  x  1(P18 

2.12  x  1(U18 

1.34  x  10  18 

Current.  Adaptive 

3.00  x  1()-'8 

3.06  x  10-18 

5.28  x  1(U18 

9.84  x  10  18 

Previous  Adaptive 

4.77  x  10-18 

5.47  x  1()-18 

3.73  x  10“ 18 

3.38  x  10  18 

Uniform 

5.38  x  10  18 

5.80  x  nr18 

3.79  x  10  18 

3.17  x  10  18 

This  table  shows  that,  in  the  lower  atmosphere,  the  current  hybrid  model  on  the  uniform 
grid  matches  the  Cf  mean  more  closely  than  Jackson’s  and  Dewan’s  model.  I11  the  upper 
atmosphere  below  20  kin,  above  which  level  grid  resolution  is  not  sufficient  to  resolve  wind 
shear,  the  current  adaptive  model  results  match  more  closely  the  observational  mean  and 
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deviation  than  the  previous  adaptive  results  and  Jackson's  model.  Examination  of  the  De¬ 
wan’s  model  results  indicates  that  the  apparent  improved  agreement  by  this  model  is  biased 
by  a  single  peak  between  11  km  <  z  <  14  km.  The  Dewan’s  model  indicates  considerably 
less  structure  than  the  adaptive  model  over  the  remainder  of  the  range.  When  improved 
grid  adaptation  scheme  is  used,  the  grid  refinement  in  the  upper  atmosphere  improves  the 
C%  prediction  over  the  previous  grid  adaptation  scheme,  which  did  not  cause  sufficient  mesh 
refinement.  But  the  mean  in  the  lower  atmosphere  is  reduced  due  to  coarser  resolution, 
caused  by  conservation  of  grid  nodes,  but  the  value  still  remains  comparable  to  Jackson’s 
model. 

5.  Conclusion 

The  first  goal  of  the  research  has  been  met  by  developing  a  proof-of-concept  NWP  code 
based  on  the  well  known  MM5  that  includes  a  new  hybrid  LES/RANS  turbulence  model 
developed  at  NCSU  for  C„  prediction  and  by  modifying  MM5  for,  and  installing,  a  version 
of  the  NCSU  3-D  dynamic  solution  adaptive  grid  algorithm  DSAGA.  The  LES/RANS  turbu¬ 
lence  model  extends  the  hybrid  NCSU  k- C  model  by  adding  two  other  model  equations;  one 
for  the  variance  of  potential  temperature,  and  the  other  for  its  dissipation  rate.  Therefore 
the  two  critical  quantities,  e  and  c#.  for  calculating  C%  can  be  directly  simulated  without 
further  simplification. 

Optical  turbulence  predictions  obtained  by  use  of  the  new  adaptive  MM5  code,  including 
the  hybrid  LES/RANS  model,  and  of  the  unmodified  MM5,  with  Dewan’s  and  Jackson’s 
models,  have  been  compared  with  two  sets  of  observations  conducted  at  the  Vandenberg  and 
Holloman  Air  Force  Bases  and  published  in  proceedings  with  sample  results  included  herein 
for  the  Holloman  case.  The  results  show  that  the  new  hybrid  model  and  the  DSAGA  adaptive 
algorithm  combine  to  predict  C„  with  variations  and  mean  values  that  approach  those  in 
the  binned  observations,  thereby  providing  more  useful  predictions  for  operational  use  as 
compared  with  information  derived  from  the  unmodified  code.  However,  this  conclusion 
must  be  qualified  by  noting  that  the  adapted  solution  is  not  yet  accurate  in  the  usual  sense 
(direct  comparison  with  the  raw  balloon  data)  as  the  shear  layers  and  other  variations  are 
displaced  and  reduced  by  boundary  effects  and  dissipation  in  the  code.  When  grid  spacing  is 
large,  it  was  noted  that  an  improved  dissipation  rate  blending  function  for  the  two  types  of 
fluctuations  should  be  explored  in  future  investigations  in  order  to  improve  the  prediction, 
especially  near  the  tropopause. 

Predictions  in  both  adaptive  and  fixed  grid  modes  are  degraded  due  to  the  fact  that 
NWP  codes,  such  as  MM5,  contain  multiple  types  of  numerical  dissipation  and,  in  general, 
use  sponge  layers  at  the  upper  boundary.  Both  of  these  features  contribute  to  smoothing 
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of  shear  layers  in  the  flow  which  then  reduces  variation  in  the  prediction.  This  effect  is 
reduced  but  not  eliminated  by  adaptation  of  the  grid.  Further  work  with  standard  NWP 
codes  should  include  efforts  to  reduce  the  amount  and  types  of  numerical  dissipation  in  these 
codes. 

The  adaptive  MM5  developed  during  this  project  is  a  research  oriented  code  only,  as 
further  development  was  halted  in  favor  of  producing  an  adaptive  version  of  WRF-ARW, 
including  the  hybrid  LES/RANS  model,  as  noted  in  the  introduction,  the  primary  code 
developer,  Dr.  Xiao,  left  the  project  just  as  work  began  on  the  WRF  version.  Due  to 
few  current  doctoral  graduates  with  CFI)  code  development  skills,  compounded  by  visa 
restrictions  for  some  candidates,  it  required  eleven  months  to  find  a  qualified  replacement  for 
Dr.  Xiao.  Because  of  this,  the  WRF  version  was  not  completed  by  project  end.  However,  the 
development  is  being  continued  under  separate  DOD  funding  and  will  therefore  be  available 
for  operational  use  when  completed. 
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Appendix:  Turbulence  Model  Equations 

A-Equation: 
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where  the  Stj  is  Kronecker  Delta,  qn  is  the  magnitude  of  the  angular  velocity  of  the  earth 
and  €ijk  is  the  permutation  tensor. 


(-Equation: 


0  d  d 

+  “  UP, 


Ht\  dC 
P  +  — 


€ mijP 


dQj 

dx^ 


d 

dxn 


(«) 


dk 

dXjyj 


pt  dQ.i  ( dQi  dQj 

ciq  J  dxj  J  oT  dxj  \  dxj  dxi 

d 5 


1  „ 


2^TiiUt-mSli  +  ^SUQjSij  +  2fcam  (y)  ; 


ku 


+  yo-Aj  +  -3jj  )  pQSij  — 

+  "  ' 


3  ”'V  13  Rk  +  S 

dk  <9(  Qj 


K 


3/2 


dxi  dxm  S 2  +  £  22/2 


+  max(Pc,  0)  -  ‘IpQSa  -  Q  0 

+  2QoC{£pr)jQj  y/c, 


Aft 

Tpk 


(up)  +c<'lS1,b'/kupup 


(A. 2) 


20 


where 
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Potential  Temperature  Variance  ( 0 "2)  Equation: 
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where 


Qt,j  -  pu"0"  =  ~p  +  Ch,g-^r-ftjzO"2j  (A. 13) 

at  =  (1  —  r)at,RANS  +  (A. 14) 

Oit.RANs  =  2  (Ch,\kro  +  Ut.RANs/ 0-89)  (A. 15) 
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Th  =  6"2lt8,kybrrd  (A.  17) 

Tg  =  0"2/ee  (A.  18) 

and  Cp  is  the  specific  heat  of  constant  pressure. 


Equation  for  the  Dissipation  Rate  of  Potential  Temperature  Variance: 
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Model  Constants 


The  constants  in  the  model  equations  are  given  below: 


Location 

Model  Constant 

Value 

k-ei  illation 

Ok 

1.8 

ck 

2.0 

C, 

0.6 

C-equation 

1.4 

Or 

0.07 

«3 

0.35 

6 

0.01 

0h 

2.37 

06 

0.10 

07 

1.50 

08 

2.3 

Q.O 

0.0 

Q.i 

1.0 

Q.2 

0.0 

Q.:» 

0.0 

1'l.RANS 

C, 

0.09 

I'tJJCS 

c. 

0.01 

r 

<*i 

10*5 

Qt.j 

ch.g 

-1.0 

(h.LES 

Ch, 

10~4 

<*t.ltANS 

Cu 

0.0648 

$"'--equat  ion 

Cf 1.2 

0.50 

chA 

-0.4 

Cu,t 

2.5  x  10“4 

ffl-equation 

Cfu  5 

-0.05 

CV« 

-0.12 

Of, 7 

1.45 

Of  ,8 

0.76 

Of, 9 

0.87 

Of, 10 

0.25 

23 


(iu)A 


-500000  0  500000 


x(m) 


Figure  1.  Domain  layout  with  terrain  elevation  contours. 
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Figure  2.  Adaptive  grid  with  weight  function  contours  and  balloon  trajectory 


25 


z  (m) 


Figure  3.  Comparison  of  potential  temperature  profiles. 
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Figure  4.  Comparison  of  mesh  spacing  profiles. 
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Figure  5.  Comparison  of  CjJ  profiles. 
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Figure  6.  Comparison  of  C £  profiles. 
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Figure  7.  Comparison  of  meridional  wind  profiles 
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